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Abstract 

We consider an oriented network in which two subgraphs (modules X and Y) , with a given intra-modular 
edge density, are inter-connected by a fixed number of edges, in both directions (nodes in X to nodes 

^—^ in y, and conversely). We study the adjacency spectrum of this network, focusing in particular on two 

Q aspects: (1) the changes in the spectrum in response to varying the intra and inter- modular edge density, 

• ^ and (2) the effects on the spectrum of perturbing the edge configuration, while keeping the densities fixed. 

O^ Since the general case is quite complex analytically, we adopted here a combination of approaches that 

allowed us to formally understand some simpler particular cases, and numerically simulate some general 

^-H results. After investigating the behavior of the mean and standard deviation of the eigenvalues, we 

^ conjectured the robustness of the adjacency spectrum to variations in edge geometry, when operating 

CN under a fixed density profile. We remark that, while this robustness increases with the size N of the 

network, it emerges at small sizes, and should not be thought of as a property that holds only in the 

large N limit. We argue that this may be helpful when studying applications on small networks of nodes 

(such as small brain macro-circuits). 

^y^ We interpret our results in the context of existing literature, which has been placing increasing attention 

^-H to random graph models, with edges connecting two given nodes with certain probabilities. Based on 

k/ similar analogies, we carry out a large N limit comparison between our network and a classical prob- 

. i-H abilistic model described in a recent paper by Narakuditi and Newman. We discuss whether / how 

^^ properties such as Wigner's semicircle law, or effects of community structure, still hold in the context 

H of a model which is density based, rather than probability based. 

Finally, we suggest possible applications of the model to understanding synaptic restructuring during 
learning algorithms, and to classifying emotional responses based on the underlying adjacency config- 
uration of the emotion-regulatory neural circuit. In this light, we argue that future directions should 
be directed towards relating such underlying hardwiring to the temporal behavior of the network as a 
dynamical system . Understanding the effects of configuration on coupled dynamics in a system in which 
each node is a nonlinear oscillator is of great importance for a wide variety of applications. 

Keywords: random matrix, oriented graph, edge density, eigenvalue spectrum, robustness. 



1 Introduction 

Eigenvalues of graphs are useful for controlling many graph properties and consequently have numerous 
algorithmic applications. Many studies have examined random graphs with a general given expected degree 
distribution, and established bounds or other descriptions/estimations of their adjacency spectra. While 
it is well know that the largest eigenvalue of such an adjacency matrix is determined by the maximum 
degree of the graph(7Ti), and the weighted average of the squares of the expected degrees (d) [^], recent 
work on random matrices has delivered more accurate estimates. For example, Chung et al. [3] have 
investigated an ensemble of random uncorrelated, undirected networks, and found that, in the large N 
limit, the expected largest eigenvalue is determined by the ratio of the second to first moment of the 
average degree distribution {(P)/{d), together with the expected largest degree dmax- More generally, for 
directed (oriented) networks without edge degree correlations, a first order approximation to the leading 
eigenvalue is given by {d^^d^^^)/{d)^ where d^^ and d^^^ are respectively the in and out-degrees of the graph, 
and (d^^) = (d^^*) = (d) [ ]. 

It is thus clear that the edge in/out degrees, as well as their correlations, have crucial effects on 
the leading eigenvalue. In general, a graph's defining feature is its distribution of edges. Among other 
properties, edge density, edge clustering, presence of hubs, have been intensely studied. Detecting and 
interpreting the modularity of a network (i.e., the presence of community structures within the graph, 
defined as densely connected groups of nodes, with sparser inter-group connections) has been recently of 
particular interest. Whether the graph represents the architecture of a web network or that of a neural 
circuit, the structure of a dynamic disease system or of a social group of friends, modularity reflects 
into adjacency properties of the network, thus modulating its structural and functional properties, and 
implicitly controlling the temporal behavior of the system. 

In this study, we focus on understanding the spectral properties of a particular type of networks. We 
consider oriented graphs with two subgraphs X and Y, each composed of N nodes. Within both X and y , 
the intra- modular edge density is fixed to a fraction 7 (out of the possible maximum of N'^). The density 
of the X-to-y edges is fixed to a fraction a of the N'^ possible X — Y connections, and the density of 
the y-to-X edges is fixed to /3 of the N'^ possible Y — X connections. Throughout the paper, we will use 
the term "modules" for X and Y (and not the more common term of "communities," since the latter has 
been used in literature in the context where the probability of two nodes within the same community to 
be connected by an edge is higher than for two arbitrary nodes which are not. In contrast, our a, /3 and 7 
take any values of the form -^ G [0, 1], where A: is a positive integer < N'^ (allowing the cases 7 < a,/3). 
In this setup, when 7 = the modules are totally disconnected, and when 7 = 1 the modules are fully 
connected (cliques). 

This graph structure may be regarded, for example, as a schematic architectural representation of 
a neural circuit. In the Discussion, we revisit an existing application, in which X and Y represent the 
excitatory, respectively inhibitory modules of a neural feedback loop, so that X projects to Y through a 
excitatory connections, and y, in turn, modulates X through /3 feedback, inhibitory connections. In such 
a circuit, the overall connectivity density may remain constant during a cognitive process such as learning, 
even though the network may exhibit high plasticity, and constantly inspect a variety of edge geometry 
combinations. Throughout the process, some connections may be lost, while others appear, as part of a 
"synaptic probing" process which has been empirically observed and is still being studied. 

The adjacency matrix of such an oriented graph is a 2N x 2N binary block matrix of the form: 

A 1 

where the blocks S and R have fixed edge density 7, while A and B have fixed edge 



B 



R 



densities a and /3, respectively. We study the spectrum of eigenvalues of this matrix, with a particular 
focus on understanding how these eigenvalues are perturbed (a) when changing the density profile (a, ^5, 7) 
and (b) when changing only the edge distribution, while keeping densities fixed. We use a combination 
of analytical and numerical methods to understand the statistics (mean and standard deviation) of this 



adjacency spectrum. 

Our work is organized as follows. In Section 2, we focus on the particular case of two interconnected 
cliques, 7 = 1 (as used in our previous modeling work). In Section 3, we generalize the connectivity scheme 
to relax the requirement of full connectedness of the modules (7 < 1), and we analyze how the robustness of 
the spectrum changes under this relaxation. In Section 4, we put our results in the context of the existing 
work on eigenspectra of random graphs. We also discuss the possible significance of our results in learning 
and neural dynamics; as an example, we revisit an application previously worked out, in which different 
brain connectivity patterns reflected into different power law behavior of the neural signals, and translated 
into different types of emotional processing. We finally propose exploring the idea of overlaying nonlinear 
dynamics on this hard-wired structure, so that the nodes of the graph are coupled oscillators. 



2 Two interconnected cliques 

We first consider the case in which the modules are fully-connected (cliques, see Figure 1). This means that 
the diagonal blocks of the adjacency matrix T are S — R — M (where M is the appropriate size matrix 
with all entries equal to one). The off-diagonal N x N blocks A and B are random binary matrices, with 
fractions a and respectively /3 of ones. By discussing the effects of edge density we mean analyzing how 
the spectrum of T changes when the values of a and /3 are varied; we will represent these changes in the 
form of surface plots with respect to pairs {a, (3) G [0, 1]^. By discussing the effects of geometry we mean 
understanding the effects on the spectrum of the edge distribution, under the constrain of fixed densities 
a and /3, respectively. We will measure these effects by estimating the mean and standard deviation of the 
eigenvalues of T over the edge geometries permissible by any fixed density pair (a,/3). 





'ys 



Figure 1: Schematic representation of the network for N = b nodes per module. Module X is shown 
on the left; module Y is shown on the right; they are both fully -connected, local sub-graphs of the full network. The 
dotted red arrows represent the X — Y connections, and the dotted blue arrows represent the Y — X connections, all 
generated randomly for low connectivity densities of a = P = .25, to maintain clarity of the illustration. 

Definition 2.1. For fixed < a, 13 < 1, we call D^'^ the distribution of 2N x 2N adjacency matrices 
T with densities with off-diagonal blocks A and B of densities a and /3, respectively. We call C^ the 
corresponding distributions of each of the eigenvalues Xj ofT (with j = 1, 2N). 

It is easy to see that the cardinality |P^'^| = C^2 C^2 • While in general the exact eigenvalues of T 
depend on the representative T E V^'^ (i.e., on the actual exact positions of the I's within the blocks A 
and S), we notice the following two boundary results (i.e., for a or /3 in {0, 1}): 

1,/3| 



Lemma 2.2. Fixing a = 1 fixes the eigenvalues ofT, so that \C- \ — 1, for all jl, 2N (i.e., the eigenvalue 



distributions are all trivial for a = 1. More precisely, the eigenvalues of any T G V^'^ are given by (from 



largest to smallest in absolute value): Ai = A^ + N\f^, \2 — N — N\f^, and A3 = . . . = \2n — (with 
multiplicity 2N — 2). Similarly, for ^ — 1, the eigenvalues of any T G V^'^ are given hy \i — N ^ N\/~a, 
\2 — N — N\/a, and A3 = (with multiplicity 2N — 2). 



Proof. We will calculate directly, for T G P^'^, the eigenvalues A and eigenvectors 
V = [vi, . . . , Vn]^ and W = [wi, . . . , Wn]^). 
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Call S^ = Zv^J ^^^ ^^ ~ /v^J' ^^^ ^3~ ^^^ i~^^ ^^^ ^^ ^^^ block matrix S, with (fj{B) being the 



number of Is in that row. We then have that: 



T^y + S^ = Xvj, for all j = 1, A^ and 



BjV + S^ = Xwj, for ah j = 1, N. 



If A 7^ 0, then vi = V2 



V, implying that T^^ — {^ ~ N)v. It follows that: ipj(B)v + S^ 



Xwj 



for all j = 1,N. By summing up, and using the fact that J2j=i ^ji^) ^ ^^/^^ we get: 

N^^v + N{X - N)v = A(A - N)v 

" T/ 1 r n 1 

. We then have that: (A - TV)^ = TV^^, hence 



Clearly i; 7^ 0, otherwise it;j = for all j, and 

A = TV ± TV V;5. 



1/ 



In conclusion: any matrix T G P^'^ has one largest eigenvalue Ai = A^ + N^/^^ with eigenvector 
Vj = V = N + TVy^, it;^- = ^j{B) + v^, and a second largest eigenvalue A2 = A^ — Ny/^, with eigen- 
vector Vj — V — N — N\/^, Wj = y^j{B) — \f^. The rest of 2N — 2 eigenvalues are zero. Note that, in the 
case of /3 = 1, then Ai = 2N and A2 = as well. □ 



Lemma 2.3. Fixing a — ^ fixes the eigenvalues of T, so that \C- \ = 1, for all j. The eigenvalues 
of any T E P^'^ are given by: Ai = A^^ A2 = . . . = A2Ar = (with multiplicity 2N — 1). Similarly, for 
/3 = 0^ the eigenvalues of any T G V^'^ are given hy Xi — N , X2 — . . . — X2N — (with multiplicity 2N — 1). 



Proof. Similar to that of Lemma 2.2. 
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When a, /3 7^ 1, the distributions C- are no longer trivial. If we restricted our interest to finding only the 



leading eigenvalue of the matrix T 
to assist us (power method, graph Lap^ 



M 
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A 



M 



, there are a variety of existing tools that we could use 



.acian). However, even the computations involved in a task such as 
expanding the powers T^ (equivalent to finding all paths of length exactly k in the graph) become very 
complex quite fast (see Appendix). It is in this light that, from this point on, we will proceed numerically 
to support a few conjectures. 

Our goal is to obtain descriptions of C^'^ for values of densities a,/3 E (0, 1); in particular, we want 
to estimate their means and standard deviations, and observe how these depend on the values of a and /3 



and on the size N of the network. For small network sizes {N < 4), the mean and standard deviation of 



the entire distribution >C^'^ can be computed directly, for each a, /3 and j (see Figure 2a and b). However, 
for larger values of N, inspecting of the whole distribution becomes computationally very expensive (e.g.. 



for TV = 5, we have |pO-48,o.48| 



(G 



12^2 

25 y 



10 configurations, although some will produce identical 



spectra). So, for larger A^s, we inspected instead only a sample y C P^'^ (we considered \y\ = 50 x 50, 
or \y\ = 100 X 100 configurations of T in P*^'^), and estimated the means and standard deviations based 
on this sample (see Figure 2b and c for a comparison between the whole-distribution and a sample-based 
results for A^ = 3). 
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Figure 2: Expected spectrum and variation of leading eigenvalue for A^ = 3; as functions of the densities 
a and f3. A. Means of the eigenvalue magnitudes, calculated formally over all T G D"'^ (i.e., over all possible 
combinatorial configurations of the blocks A and B), for each given pair of densities (a,/3). B. Standard deviations 
of the eigenvalue magnitudes, also calculated formally over all possible combinatorial configurations. C. Standard 
deviations of the eigenvalue magnitudes, calculated using a sample of 50 representative matrices for A, and 50 for 
B, for each given pair of densities (a,/3). 



2.1 Estimating the large eigenvalue means 

There are a few contexts in the literature on eigenspectra of random graphs that relate to our problem, 

among which a known result due to Juhasz [ ]. Viewed in the reference's general framework, the adjacency 

1 a 
matrix T is a block matrix with (weighted) density matrix D 



/? 1 



, whose eigenvalues are /xi^2 



1 ± ^Ja^. According to the main theorem in the referenced paper, T has 2 eigenvalues that are large in 
absolute value (of order A^), with: Ai^2 = N ^N \fa^^o{N^'^^^') in probability, while the other eigenvalues 
are of order o{N^I^^^^ in probability (for any 6 > 0). 

Based on this result, \\^2 = A^ ± N \Jol^ would be the first candidates to check, as possible expressions 
for the means of the two large eigenvalues of T in terms of the densities a and ^. The formula looks 
particularly promising, since it seems to naturally extend the boundary expressions obtained in the two 
lemmas for a or /3 = 1 in {0, 1}, and since it matches almost exactly the surfaces we obtained numerically 
for Ai^2- However, this was quickly refuted by the direct computation of the mean and standard deviation 
for small network sizes (A^ = 3,4), although the formula seems to give a very close estimate, which might 
be exact in the limit as A^ ^ oc. This leaves us with the question of computing and/or understanding the 
source of the small correction terms. 



2.2 Spectrum variance 

For a fixed A^, we noticed that, as expected, the eigenvalue distributions get thicker at intermediate values 
of a and /3. Figures 2 illustrates the surface plot (with respect to a and /3) of the standard deviation of 
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Figure 3: Expected spectrum and variation of eigenvalues for N = 7, estimated numerically by considering a 
sample of 50 matrices A and 50 matrices B, for each pair of densities (a, 13). A. Surface plots of expected eigenvalue 
magnitudes. There are two leading eigenvalues, whose surfaces are very close to the surfaces N ± N^/ap; the other 
eigenvalues are close to zero. B. Surface plot of the standard deviation estimate for the largest eigenvalue, for each 
pair of densities (a, /3). B. Surface plot of the standard deviation estimate of the second largest eigenvalue magnitude. 
B. Surface plot of the standard deviation estimate of the third largest eigenvalue magnitude. 



the leading eigenvalue of T, computed over the whole distribution (Figure 2b), as well as over a random 
sample of 50 x 50 representatives for each pair (a,/5). Similarly, or Figure 3, we used N — 1 and a lattice 
of 50 X 50 matrices for each pair (a, /3) to estimate numerically the standard deviations of the highest three 
eigenvalue magnitudes. In both cases, the surfaces decrease from the center towards the edges, illustrating 
the narrowing of the eigenvalue distributions when a or /3 get closer to and 1. 

Following the same logic ("higher cardinality likely produces higher variance"), one would have expected 



iQ;,/3| 






the eigenvalue standard deviations to increase when the size N is increased (recall that \D' 
which increases greatly with N). Juhasz' estimate seems in tune with this, allowing a correction term of 
magnitude o(A^^/^+^). In reality, the "widths" of >Cj'^ do not seem to increase as \/]V, but rather saturate 
as N increases (the case of the small eigenvalues) , or even narrow down (the case of the leading eigenvalue 
distributions). In Figures 4a,b,c we illustrate, as surface plots, the behavior of the standard deviations for 
the three leading eigenvalues, for 2 < A^ < 7; in Figure 4d we show that the maximum standard deviation 
(as approximated by our algorithm) appears to be decreasing as a function of size N . Based on our plots, 
we suspect that, for any pair (q;,/3), the standard deviations change slower with N SiS N increases. Com- 
puting these rates, however, seems quite difficult, since the distances between two consecutive approximate 
surfaces are hard to estimate numerically with enough accuracy, even for small values of N. 
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Figure 4: Illustration of the evolution of the standard deviation of the eigenvalue magnitudes, when 
increasing the network size N, A. Each surface plot represents the standard deviation, with respect to {a^P), 
of the leading eigenvalue of T, for a different size N , from N = 2 (lightest surface, top) to N = 7 (darkest surface, 
bottom). B. Each surface plot represents the standard deviation of the second leading eigenvalue magnitude, for a 
different size N, from N = 2 to N = 7 (notice that the surfaces cross each other). B. Each surface plot represents, 
for a different size 2 < N < 7, the standard deviation for a small eigenvalue of T (the third in magnitude, in this 
case). For convenience, we plot here one minus the actual value; the figure is vertically flipped, in order to make all 
surfaces visible. D. Illustration of how the standard deviations of the eigenvalue magnitudes evolve when the size N 
is increased. We plot the maximum value attained on the surface corresponding to each N , for a few representatives 
of the spectrum (for the too eigenvalues large in magnitude, as well as for the third and the last in magnitude). 



We summarize our theoretical and numerical observations in this case in the form of a conjecture, which 
remains open to a more rigorous investigation: 

Conjecture 2.4. In the case of fully- connected modules 7=1 (i.e., 5 = i? = M), the spectrum of the 
matrix T varies with respect to the inter-modular densities a and /3 of the blocks A and B as follows: 



(i) The spectrum has two eigenvalues Ai and A2 whose means are large in magnitude, while the other 
2N — 2 have means which are small in magnitude (close to zero). 

(a) For each size N and each density pair {a^ji), the mean magnitudes of the two leading eigenvalues 
(over all adjacency configurations corresponding to (a, P)), are given approximately by N ± N^fap. 



(in) For any size N, the standard deviation of each eigenvalue magnitude is a ''unimodaU^ surface, with a 
point of maximum in the open square (0, 1)^^ and which is zero when a = 0, 1 or /3 = 0, 1. 

(iv) As N increases, the surface- maxima saturate asymptotically. The change with N (difference between 
two consecutive maxima) decreases as N gets large. 

(v) For the leading eigenvalue Xi, the standard deviation maximum is very small (one degree of magnitude 
smaller than the mean Xi, or than the corresponding maxima for the other eigenvalues). Moreover, 
the standard deviation of Xi decreases monotonically with N, for each fixed pair (q;,/3). 

3 Two interconnected modules 

In this section, we relax the "fuh-connectedness" condition for the two modules, to explore other intra- 
modular edge configurations. 
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Figure 5: Illustrations of the standard deviations of the eigenvalue magnitudes for different values of intre-modular 
connectivity j, for N = 3. A. Each surface represents the standard deviation of the leading eigenvalue, for one value 
ofje [0, 1] (low to high surfaces, as 7 decreases). B. Cross-section of the surfaces shown in {A), for fixed a = 5/9. 
C. Similar cross-sections, shown for the second leading eigenvalue. D. Similar cross-sections, shown for a small 
(third in magnitude) eigenvalue. 

In Figure 5, we illustrate the dependence of the standard deviations not only on the inter-modular 
edge densities a and /3 (represented on the x and ?/-axes), but also on the intra-modular density 7 (in each 



panel, different plots correspond to different values of 7 G [0, 1]). In Figure 5a, we represented the standard 
deviations of the leading eigenvalue, for A^ = 3 and 7 < 1, as surface plots with respect to pairs («,/3). 
(For better visualization of the surfaces, we omitted the boundaries (a, 1) and (l,/3), where the standard 
deviations are zero). Notice that, for pairs (a,/3) close to the corner (1, 1), the surfaces are (as one would 
expect) barely affected by the value of 7, as long as 7 > 0. More generally: although the stdev surfaces 
do not generally exhibit the same shape and unique centrally-placed maximum as in the particular case of 
Section 2, the unimodality is preserved in cross-section (in Figures 5b, c and d we show such cross-sections 
for the first three largest eigenvalues, for A^ = 3). 

As 7 decreases, the standard deviation surfaces raise higher. However, a main property appears to be 
preserved from the fully-connected case 7 = 1 to the more general case of 7 > 0: the standard deviations 
don't escape to oc as A^ increases, but are rather bounded by a fairly small upper bound (close to 1). In 
the case 7 = 0, however, this property breaks down. It is interesting to notice that the standard deviation 
for the leading eigenvalue, which is very small when 7 = 1, increases slowly as 7 decreases from 1, then 
faster as the values of 7 get close to 0; the complete crash, however, only happens at 7 = (also see Figure 
6). In other words, the full-connectedness of the two moduli confers very high spectral robustness to the 
network, which is surprisingly well preserved as one loosens this condition, pruning out random edges and 
lowering the intra-modular density. The property is completely lost, however, when there are none of the 
intra-modular edges left to connect the nodes. We further interpret this in the Discussion. 
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Figure 6: Behavior of eigenvalues of T a^ 7 decreases from 1 to 0. Here, N = 3, and a = 5/9. In each panel: 
7 = 1 (panel A), 7 = 4/9, (panel 'B), 7 = 2/9 (panel C) and 7 = (panel T>). The curves represent the expected 
eigenvalue magnitudes, and the error bars represent the respective standard deviations. 



If the graph is the underlying hardware for a dynamical system, the local linearization (Jacobian) matrix 
around the system's equilibria is likely to depend, among other parameters, on the adjacency spectrum, in 
particular on the separation of eigenvalues. Here, we also illustrate how the eigenvalues are separated, and 
how this separation if affected when using the different connectivity schemes we considered throughout 
this section. Figure 6 shows the means of the eigenvalue magnitudes together with the respective standard 
deviations, for four level of intra-modular connectivity (each panel corresponds to a different value of 7, 
as described below); we fixed a = 5/9 and allowed /3 to vary in each case along the x-axis. When 7 = 1 
(Figure 6a), the standard deviations are small (as shown in Section 2), and the first and second eigenvalues 
(whose means are well approximated in magnitude by N ^N^fa^) remain to a large extent separated from 
the other small eigenvalues. As 7 decreases from 1, this situation gradually changes, and the large expected 
eigenvalues decay as A^7 ± Ny/a^ (Figure 6b, c), to eventually completely collapse at 7 = (Figure 6d). 

4 Discussion 

4.1 General comments 

In this study, we have investigated, using analytical and numerical computations, the adjacency spectrum 
of an oriented graph, in which the nodes of two modules connect though fixed numbers of randomly placed 
edges within each module, as well as across modules. We concluded that, when fixing the number of both 
intra and inter-modular edges, the adjacency spectrum of the network remains in general sufficiently robust 
under particular edge configurations (geometries), suggesting that simple algorithms in such a system may 
also remain unaffected by constrained geometry changes. 

There is a very large body of work addressing properties of random matrices [23], whose entries are 
drawn independently out of a given (typically normal) probability distribution. If, in addition, the matrix 
represents the adjacency of a random graph, so that each entry equals 1 with a given probability, there 
are classical methods used when looking for properties of the the spectrum, such as spectral radius, or 
spectral density. Our model differs from most of these approaches in that it conserves the number of edges 
within/between modules, rather than fixing the probability of independently having an edge that connects 
two given nodes in the same/different modules. In our setup, the entries of the adjacency matrix are neither 
independent, nor identically distributed. However, while classical results (such as Wigner's semicircle law) 
require the entries to be identically distributed, various extensions have been worked out, for models which 
don't necessarily abide by these properties. 

Consider, for example, the configuration model [ ], whose spectral properties have been addressed by 
numerous studies. Since its edges are not statistically independent, a direct analytical approach is very 
difficult; existing results range from approximating the full spectrum [u], to formally deriving the expected 
values of the leading eigenvalue, but only in the large N limit ["]. In a recent paper, Newman et al. [ ] took 
an indirect approach: they considered instead a model with the same degree sequence as the configuration 
model, but in which the number of edges between any two nodes was drawn independently from a Poisson 
distribution. They then showed that the spectra of the two models agree in the large N limit. 

Below, we illustrate the same idea, by carrying out a large N limit comparison between our model and 
a corresponding model, with independent, stochastic edges, considered by Nadakuditi and Newman [ 4]. 
In the referenced paper, they considered a stochastic, non-oriented network with two communities, and 
computed the ensemble-means for the two large eigenvalues of its symmetric adjacency matrix, in the 
large N limit. The method involved first finding the eigenvalues of the modularity matrix, then showing 
that these are identical in the large N limit to the eigenvalues of the adjacency matrix. Their asymptotic 
expressions zi and Z2 were computed in terms of cin = npi^ and Cout — ^Pout (where, with the notations in 
the original text, n is the matrix size, pi^ is the probability of two nodes within a module to be directly 
connected, and pout is the probability of two nodes which are not in the same module to be directly 
connected). More precisely: 
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Figure 7: Comparison between our results and those of Narakuditi- Newman [14], in the case of a two 
modular, non-oriented graph. We compare the values of zi^2 (solid corves in purple and brown, respectively) 
with the formal means Ai^2 (dotted curves in blue and green), and their close approximations 7V± — (solid curves in 
yellow and cyan). A. Comparison for N = 4. B. Comparison for N = 100; here, we used only the approximations 
N ± —, since the formal means are computationally extremely expensive. 
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Z2 = 7:(<^in - Coutj + 






With our notation, Cin = 2N, Cout — "^a/N and the adjacency matrix is symmetric (/3 = a). Accounting 
for the absence of loops in the referenced model, which are present in ours, we get: 



zi 

Z2 



N^ + a 

N 
N'^ -a N^ ^a 



N 



N^- 



1 



a 



so that zi > Z2 ii oi < N{N — 1). In Figure 7, we show a comparison between our results and those 
of Nadakuditi and Newman, when applied to two fully-connected communities, by illustrating on the 
same axes zi^2^ the formal means Ai^2 5 cind their close approximations obtained earlier as A^ it — -. The 

approximations get exact in the large N limit, at least for values oi a < N'^ — N\/2N ^ where zi has its 
global minimum (after which it shoots up, detaching from the graph of A2). 

Lastly, the reference remarks on spectral distribution of the modularity matrix (i.e., in the large N 
limit, also the spectral radius for the adjacency matrix of the non-oriented graph studied). The spectrum 
consisted of a continuous semicircular band of eigenvalues, plus the singular leading eigenvalue. As long 
as the leading eigenvalue is well separated from the semicircular band, there is evidence of community 
structure in the network; when the leading eigenvalue passes the edge of the band [z — \/c\^ + Cout)? the 
community structure is no longer detectable. As we have already suggested in Section 3, the property 
extends to the case of the oriented random graph that constitutes our study case. In Figure 8 we show, 
for A^ = 4, a sample (100 configurations) of the spectrum, observing the separation between eigenvalues, as 
the modularity of the network changes. Since the adjacency matrix is no longer symmetric, the eigenvalues 
are plotted in the complex plane. All eigenvalues are distributed (uniformly?) within the unit disc, except 
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the first two largest in absolute value, which, for 7 = 1, are real and significantly larger than 1. When the 
"community structure" decreases (e.g., 7 decreases), first the second largest eigenvalue collides into the 
unit disc, and starts diffusing around its boundary. If we continue decreasing 7, the leading eigenvalue will 
also become indistinguishable from the pool distribution. 



iV = 4, a = 0.25, /3 = 0.75 




iV = 4, a = 0.5, /3 = 0.5 




- ^ 




0^ 





Figure 8: Illustration of separation of eigenvalues, when changing the community structure. The 

eigenvalues are plotted in the complex plane: the leading eigenvalue in absolute value is shown in red, the second 
largest in green, the rest in blue. All plots are for N = 4, and are based on samples of 100 matrix configurations, 
under the following restrictions: A. a = 1/4, j3 = 3/4^ 7 = 1,3/4,1/2,1/4 (from top to bottom). B. a = 1/2^ 
/3 = 1/2, 7 = 1, 3/4, 1/2, 1/4 (from top to bottom). 



4.2 Applications to learning 

In summary: while a lot is known about certain classes of random matrices (e.g., with normal, or in- 
dependent entries), little is known about other types of ensembles. Normal distributions are ubiquitous 
in nature; however, natural processes may involve many other sources of noise and of randomness, and 
studying them may become equally useful. 

For example, let's imagine the oriented graph considered in this paper to represent a network of coupled 
neurons, so that each edge represents a synapse with a corresponding synaptic strength, or "weight " (so 
that the "effective connectivity matrix" of such a network would be a weighted mutation of the adjacency 
matrix). Synaptic updating has been known to be for a long time the physiological basis of learning, but 
the exact ways in which such a process is implemented biophysically are still under investigation. 

There are many different models describing, qualitatively or quantitatively, the synaptic adjustments 
that may take place in such a network, during a process such as learning. In general, the process is 
suspected to involve not only weight changes of existing synapses, but also activation of "silent" sites. 
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creating of new connections, or pruning of the existing ones. In terms of our model, this means that not 
only the edge weights, but also the edge distribution might change during learning. 

A clear biological restriction to this updating game has to be that the synaptic strengths are somehow 
prevented to increase without bound, which is why most models incorporate a normalization scheme. 
However, again, there is very little agreement on how this normalization step is actually implemented by 
the brain. Some rules assume the process to be local (subtractive normalization [15, 9, 26, 12], a weight- 
dependent rule [i ] or a BCM rule [ ]), but one can imagine various other ways of insuring stability, possibly 
involving "homeostasis" or "synaptic scaling" [^^4, '^^. Many models support a global normalization, for 
which the state of the whole network is assessed at each updating step, and a specific norm is imposed at 
each weight update. 

In this light, it becomes important to understand possible normalization mechanisms for optimal weight 
updating. While most models of learning introduce the updates into the weights themselves, the brain 
may normalize (at least in the short term) by simply maintaining the overall number of active connections 
in the network approximately constant, so that, in the updating process, loosely one synapse will turn off 
whenever a new site is activated. One would then want to understand how this architectural dynamic may 
promote/influence learning, and how the effects of geometry updating complement or compare with the 
effects of direct weight updating. Let us finally note that adding or deleting edges based on a probabilistic 
process (as described in Section 4.1) would be a local mechanism, equivalent to our model in the large N 
limit. However, even for relatively large values of N (see Figure 7b) for N = 100), the differences in the 
spectra are quite substantial. Since many brain networks plausibly operate with a few thousand neurons, 
it is safe to understand the apparent distinctions between the two models, for finite values of N. 

Knowledge of the geometry of the network is very important when determining which connectivity 
schemes are plausible to use for models of learning. Choices currently used in modeling range from 
considering fully- connected to fully-disconnected interacting modules, or layers [ "]. In our framework, 
convergence (learning) is not a priori prevented in either case. In developing future iterations of this 
model, it will also be important to explore how the learning process itself shapes the connectivity scheme. 
Siri et al. [ ] suggest that the structure emerging during learning breaks down into different numbers of 
hub-like subnetworks; this might affect the spectral robustness demonstrated in our modular. Understand- 
ing the source and limits of this robustness may be an instrument that could be used to investigate which 
architectures favor convergence under particular learning algorithms, and which not. 

4.3 Random coupled oscillators and applications to brain connectivity 

Recent studies have used graph theoretical approaches not only to understand learning, but also to inves- 
tigate more general organizational principles of brain networks [ ^, 22, 20]. With nodes and edges defined 
according to modality appropriate scales [ ], there studies support certain generic topological properties 
of the human brain architecture, such as modularity, small- worldness, the existence of hubs and other 
connectivity density patterns [10]. These properties, if proven consistent with certain physiological, be- 
havioral or genetic factors, may be usable to better understand neural processes, or even as biomarkers for 
behavioral traits or neuropsychiatric conditions. 

In previous work, we have used a graph-theoretical framework to study how network density can affect 
the complexity of signal outputs in a control meso-circuit (the prefrontal-limbic system, viewed as a network 
of hemodynamic nodes). As in our theoretical framework, we used two interacting modules (one excitatory 
- amygdala, and one inhibitory - prefrontal cortex), thus representing schematically the prefrontal-limbic 
system. While conceived to be sufficiently general to maintain relevance for a variety of classes of neural 
control circuits, this setup wa appropriate for informing our imaging results on this system in humans, 
and draw conclusions on how prefrontal-limbic connectivity may drive arousal dynamics and emotional 
responses. Woth each of the N nodes in either module acting as an stochastic linear damped (Ornstien- 
Uhlenbeck) oscillator, we studied how connectivity complexity can determine and modulate scale-free 
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features of the output signals. When ahowing the densities a, /3 to vary within the interval (0,1], the 
frequency profiles of the solutions (measured in our case as the slopes of the log-log power spectra) shifted 
closer of further from white to pink to brown noise, in a consistent way which allowed us to suggest a testable 
framework for interpreting similar power spectral trends found in our fMRI time series. Individuals with 
average emotional reactivity represent well-regulated control systems, in which excitatory and inhibitory 
influences are balanced (these individuals exhibited fMRI signals close to pink noise in both amygdala 
and prefrontal regions). Anxious individuals have relatively weaker inhibitory feedback inputs from the 
prefrontal cortex (primarily driving amygdala signals closer to white noise). Less reactive individuals have 
relatively stronger excitatory inputs from the amygdala, producing stronger feedback (inducing more white 
noise primarily for the prefrontal cortex). 

Such modeling results, in agreement with empirical finds [ ] and in with existing theoretical work [1], 
support the necessity of two key conditions for optimal function in brain networks: first, the requirement 
for a well-balanced adjacency matrix. The corresponding bidirectional graph should appropriately combine 
robust features (e.g., fully connected populations) and random edges (e.g., due to synaptic probing) so as to 
allow some flexibility, yet also render sufficient stability for convergence during a cognitive process such as 
learning. Second, there should be well-balanced connection strengths (i.e., the weights on each adjacency 
edge), hence an efficient connectivity matrix driving optimal dynamics in the system. 

In future work, we propose extending such analyses to understand interactions between nonlinear nodes. 
While biologically more plausible, nonlinearity also produces richer temporal trajectories. The dynamics of 
many types of coupled nonlinear oscillators have been widely investigated, in particular when the tempo- 
ral architecture of the system is overlayed on a hard-wired connectivity scheme between nodes/ variables. 
It has become clear that even trivial interaction schemes, in conjunction with nonlinear behavior, may 
produce highly complex phenomena. For example, symmetric weak coupling between two Wilson-Cowan 
excitatory/inhibitory units has been shown to produce interesting phase-space transitions (bifurcations) 
between stable equilibria, cycles and invariant tori, when varying the system parameters (e.g., the coupling 
strengths). In our framework, the aim would be not necessarily to locate these transitions in the param- 
eter space for each adjacency configuration, but rather to understand how the bifurcations are affected 
statistically by changes in the edge configuration. 
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